{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1d72bf50",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using the Community license in this session. If you have a full Xpress license, first set the XPAUTH_PATH environment variable to the full path to your license file, xpauth.xpr, and then restart Python. If you want to use the FICO Community license and no longer want to see this message, set the XPAUTH_PATH environment variable to: /home/maged/anaconda3/envs/balance/lib/python3.8/site-packages/xpress/license/community-xpauth.xpr\n",
      "NB: setting XPAUTH_PATH will also affect any other Xpress products installed on your system.\n"
     ]
    }
   ],
   "source": [
    "try:\n",
    "    import hsbalance as hs;\n",
    "except ImportError:\n",
    "    !pip install hsbalance\n",
    "    import hsbalance as hs;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f622ade9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1045eff9",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "\n",
    "In these notes, we are going to discuss the issues related to robust optimization.\n",
    "a. What is Robust Optimization?\n",
    "Optimization process is not determinant, this means that we have no clue or determinant check if our solution is the optimum one. Optimization is prone to error if outliers exist. Noisy data and uncertainty readings. Methods, especially least squares method, can be more effected with outliers than others.  \n",
    "b. Rotor Balancing Robust Optimization:\n",
    "In our problem of rotor balancing, solution can be greatly effected by outliers for the following reasons:  \n",
    "1. Measurements are low in numbers: Maximum field measuring points number is still tiny mathematically-wise.\n",
    "2. Sensor faults: a sensor can go wrong easily for a turbine that runs for 20 years. Calibration is not an easy or a cheap process.  \n",
    "3. Accumulation of errors: The process of measuring a point in a turbine/compressor usually involves a sensor and long cables in series (vibration cards and junctions). Those produce errors in series and cause the summation to infiltrate.\n",
    "4. Measuring times: Turbo-Machinery is often connected to large process, which make it start and stop time costly. This encourages for accurate balancing with the fewest number of starts."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab28df86",
   "metadata": {},
   "source": [
    "# Simulating Data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46553bd3",
   "metadata": {},
   "source": [
    "We generate random data with known solution W"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "dde2a986",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "real = np.random.rand(5,3)\n",
    "imag = np.random.rand(5,3)\n",
    "alpha = real + imag * 1j\n",
    "# Assume W is equal to 10 grams @ 0 degree at all planes\n",
    "W = [[10], [10], [10]]\n",
    "# Create intial condtioning column array A with random noise to simulate actual problem\n",
    "A = - alpha @ W + np.random.normal(1, scale=0.1, size=(5,1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9ebf3c96",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9546"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Error due to generated noise\n",
    "hs.rmse(alpha @ W + A)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52cf1293",
   "metadata": {},
   "source": [
    "# Create a hsbalance model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "92a68852",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([['9.448 @ 4.1'],\n",
       "       ['9.49 @ 0.1'],\n",
       "       ['9.842 @ 1.5']], dtype='<U11')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_alpha = hs.Alpha()\n",
    "model_alpha.add(alpha)\n",
    "model = hs.LeastSquares(A, model_alpha)\n",
    "W_model = model.solve()\n",
    "hs.convert_cart_math(W_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02ac52e5",
   "metadata": {},
   "source": [
    "Calculate the Error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ffe9a0a8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5666"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hs.rmse(W_model - W)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd03ba24",
   "metadata": {},
   "source": [
    "# Faulty Sensor"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa10eb68",
   "metadata": {},
   "source": [
    "To simulate a faulty sensor we will assume that the sensor at the first plane has gone completely fault and gives no output (i.e. equal 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b716ba39",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([['11.739 @ 18.8'],\n",
       "       ['5.966 @ 294.0'],\n",
       "       ['7.123 @ 35.6']], dtype='<U13')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[0] =0\n",
    "model_alpha = hs.Alpha()\n",
    "model_alpha.add(alpha)\n",
    "model = hs.LeastSquares(A, model_alpha)\n",
    "W_model = model.solve()\n",
    "hs.convert_cart_math(W_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d37dbcbf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.3969"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hs.rmse(W_model - W)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d54e25e0",
   "metadata": {},
   "source": [
    "As we see, the error has gone extremly high and obviously unacceptable.  \n",
    "If we acknowldge the sensor error, we can use weight least squares method to solve the problem instead.  \n",
    "This can be done by assign the argument solver='WLS' in solver method of least squares. and pass the coefficient column array C to the model.   \n",
    "In our case we will give all the sensor the same weight as 1. Only the sensor suspected we will assume a weight of 0. The previous setup will completely discard the sensor from caluclations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "47055e19",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([['9.52 @ 4.1'],\n",
       "       ['9.317 @ 0.4'],\n",
       "       ['9.865 @ 2.1']], dtype='<U11')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C = np.ones(A.shape)\n",
    "C[0] = 0\n",
    "model = hs.LeastSquares(A=A, alpha=model_alpha, C=C)\n",
    "W_model = model.solve(solver='WLS')\n",
    "hs.convert_cart_math(W_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "419019c2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6433"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hs.rmse(W_model - W)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa833954",
   "metadata": {},
   "source": [
    "We can notice that we restored a low level of error when discarding the faulty sensor."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c60df44d",
   "metadata": {},
   "source": [
    "# Huber Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea02b46c",
   "metadata": {},
   "source": [
    "The previous example assumes we know for sure the defected sensor, so we down-weighted it to decrease the leaking noise it causes. A more convenient method but more conservative. The method make the model less sensitive to outliers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ff0826c8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([['10.859 @ 19.6'],\n",
       "       ['7.188 @ 316.5'],\n",
       "       ['6.926 @ 16.3']], dtype='<U13')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = hs.LeastSquares(A, model_alpha)\n",
    "W_model = model.solve(solver='huber')\n",
    "hs.convert_cart_math(W_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "9bd3c129",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.8004"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hs.rmse(W_model - W)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "264ad4b9",
   "metadata": {},
   "source": [
    "The error improved but in a more conservative way."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67239fd1",
   "metadata": {},
   "source": [
    "| Method | RMSE |\n",
    "|--|--|\n",
    "|Least_squares|0.5666|\n",
    "|Least_squares with errorous sensor|6.3969|\n",
    "|Weighted_Least_squares|0.6433|\n",
    "|Huber_Robust|4.8004|"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08ab00a7",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:balance] *",
   "language": "python",
   "name": "conda-env-balance-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
